# Analysis conducted using R version 4.0.2 on OS 10.14.6
rm(list=ls())
setwd("~/Dropbox/Writing/Military Officers and Public Opinion/Data Analysis")
data <-read.csv("DataOfficers_Replication.csv", fileEncoding = "UTF-8-BOM")

install.packages("doBy")
install.packages("gmodels")
install.packages("ggplot2")
install.packages("xtable")
install.packages("XLConnect")
install.packages("xlsx")
install.packages("aod")
install.packages("MASS")
install.packages("mvtnorm")
install.packages("xtable")
install.packages("stargazer")
install.packages("dplyr")
install.packages("descr")
install.packages("DescTools")
install.packages("ggpubr")
install.packages("purrr")


library(doBy) 
library(plyr)
library(dplyr)
library(gmodels) 
library(ggplot2)
library(xtable)
library(XLConnect)
library(xlsx)
library(sandwich)
library(lmtest)
library(aod)
library(MASS)
library(xtable)
library(stargazer)
library(descr)
library(DescTools)
library(ggpubr)
library(scales)
library(rlang)

#----------------------
#Variable Creation
##Create Rank Variable

data$Rank <- NA
data$Rank[data$O1==1] <- 1 #n=3
data$Rank[data$O2==1] <- 2 #n=11
data$Rank[data$O3==1] <- 3 #n=65
data$Rank[data$O4==1] <- 4 #n=49
data$Rank[data$O5==1] <- 5 #n=36
data$Rank[data$O6==1] <- 6 #n=20
data$Rank[data$O7==1] <- 7 #n=1

FieldGradePlus <- (sum(data$Rank == 4)+sum(data$Rank == 5)+sum(data$Rank == 6)+sum(data$Rank == 7))/185 #Presented in main text as percentage of officers in rank of Maj/LCDR or above

#Create Air Force Variable (Used for Chi-Square tests; see Appendix 2.1)

data$AirForce <- NA
data$AirForce[data$Mil_Branch==5] <- 1 #n=103 Air Force, n=82 not Air Force
data$AirForce[data$Mil_Branch==4] <- 0
data$AirForce[data$Mil_Branch==3] <- 0
data$AirForce[data$Mil_Branch==2] <- 0
data$AirForce[data$Mil_Branch==1] <- 0

#----------------------

#Main Findings

##Recommendations Toward Use of Force by Condition (Findings presented in first paragraph of findings section.)

RecommendationTable <- table(data$PublicOpinion, data$Recommendation)
PropRecommendationTable <- prop.table(RecommendationTable,1)
rownames(RecommendationTable) <-c("Public Supports Military Action", "Public Opposes Military Action")
colnames(RecommendationTable) <-c("Recommend Operation", "Recommend Against Operation")
prop.test(RecommendationTable)
chisq.test(RecommendationTable) #X-squared = 9.1659, df = 1, p-value = 0.002466

##Stacked Plot: Recommendations Toward Use of Force by Condition (Not included in text/appendix)
RecommendationPercentage <- apply(RecommendationTable, 1, function(x){x*100/sum(x,na.rm=T)})
SupportPlot <- barplot(RecommendationPercentage, xlab="Experimental Condition", ylab="Percent")

##Regression Analysis (Appendix, Table 1)
m1_Recommendation <- lm(Recommendation ~  PublicSupports, data=data) 
summary(m1_Recommendation)

m2_Recommendation <- lm(Recommendation ~  PublicSupports + Gender + Age + Mil_Branch + Combat + Status + Rank + Staff + Education + Race + Political_ID + Income, data=data) 
summary(m2_Recommendation)
stargazer (m1_Recommendation, m2_Recommendation)

#Branch Effects and Chi Square Test (Appendix, Table 2)
AirForceTable <- table(data$AirForce, data$Recommendation)
PropAirForceTable <- prop.table(AirForceTable,1)
rownames(AirForceTable) <-c("Not Air Force", "Air Force")
colnames(AirForceTable) <-c("Recommend Operation", "Recommend Against Operation")
prop.test(AirForceTable)
chisq.test(AirForceTable) #X-squared = 0.0014601, df = 1, p-value = 0.9695

#-----------------------
#Justifications Analysis (Main Text, Table 1)

#Justifications: "Write a sentence or two telling us why you choose this recommendation"?
## 1- Recommendation Accounts for Public Opinion; 2- Public opinion is beyond purview of military officers; 3- Defends American interests and/or reputation; 4-Risks of launching retaliation is high; 5- No Response/Need More info to answer

JustificationTable <- table(data$PublicOpinion, data$JustificationID)
PropJustificationTable <- prop.table(JustificationTable, 1) 
print(PropJustificationTable)*100

## Intercoder Reliability Check
sum(data$JustificationID == data$JustificationID_2)/length(data$JustificationID) #0.9135135

#----------------------
# How much do military officers consider public opinion (Percentages presented in main text)

##"To what extent do you consider public opinion when making decisions on use of force?"
##(Coded on five point Likert scale from "Not at all (1)" to "Very Large Extent (5)")

ConsiderPublicOpinionTable <-table(data$ConsiderPublicOpinion)
prop.table(ConsiderPublicOpinionTable) 

# Comparison of mean levels of consideration to ensure no difference in levels of consideration in public support across experimental conditions 
##Full sample
mean(data$ConsiderPublicOpinion) #2.664865
## Means, sd, and se by treatment
Mean_ConsiderPublicOpinion <- ddply(data, c("PublicOpinion"), summarise, mean = mean(ConsiderPublicOpinion), sd = sd(ConsiderPublicOpinion), sem = sd(ConsiderPublicOpinion)/sqrt(length(ConsiderPublicOpinion)))
print(Mean_ConsiderPublicOpinion) #2.624 (.085), 2.70 (.084)

#----------------------
# Reasons for considering public opinion (Main Text, Table 2 and Appendix, Table 6)

# How important are each of the following reasons for considering public opinion?

# Main Text, Table 2/Appendix Section 2.3 and Appendix Table 6
## Note: In table 2, I combine "not at all important" and "slightly important" into a single column and "very important" and "extremely important" into another. The full table is in Appendix, Section 2.3. 

## "Public support is needed to carry out successful military operations"
OpSuccessTable <-table(data$OpSuccess)
prop.table(OpSuccessTable)

Mean_OpSuccess <- ddply(data, c("PublicOpinion"), summarise, mean = mean(OpSuccess), sd = sd(OpSuccess), sem = sd(OpSuccess)/sqrt(length(OpSuccess)))
print(Mean_OpSuccess) #Public Supports: 2.694118; Public Opposes: 2.780000)

Mean_OpSuccess_Full <- (mean = mean(data$OpSuccess)) #2.740541
sd_OpSuccess_Full = sd(data$OpSuccess) #1.072244
se_OpSuccess_Full = sd(data$OpSuccess)/sqrt(length(data$OpSuccess)) #0.07883291

m1_OpSuccess <- lm(OpSuccess ~ PublicOpinion, data=data) 
summary(m1_OpSuccess) # p-value: .5886 

## "Taking actions that run counter to public opinion jeopardizes the military's reputation"
MilReputationTable <-table(data$MilReputation)
prop.table(MilReputationTable)

Mean_MilReputation <- ddply(data, c("PublicOpinion"), summarise, mean = mean(MilReputation), sd = sd(MilReputation), sem = sd(MilReputation)/sqrt(length(MilReputation)))
print(Mean_MilReputation) #Public Supports: 2.8; Public Opposes: 2.8

Mean_MilReputation_Full <- (mean = mean(data$MilReputation)) #2.8
sd_MilReputation_Full = sd(data$MilReputation) #1.015103
se_MilReputation_Full = sd(data$MilReputation)/sqrt(length(data$MilReputation)) #0.07463188

m1_MilReputation <- lm(MilReputation ~ PublicOpinion, data=data)
summary(m1_MilReputation) # p-value:  1

## "The military carries out the policies of elected leaders who are responsive to public opinion. The military should anticipate/respond to these leaders"
LeadersPreferenceTable <-table(data$LeadersPreference)
prop.table(LeadersPreferenceTable)

Mean_LeadersPreference <- ddply(data, c("PublicOpinion"), summarise, mean = mean(LeadersPreference), sd = sd(LeadersPreference), sem = sd(LeadersPreference )/sqrt(length(LeadersPreference)))
print(Mean_LeadersPreference) #Public Supports: 3.63; Public Opposes: 3.69

Mean_LeadersPreference_Full <- (mean = mean(data$LeadersPreference)) #3.659459
sd_LeadersPreference_Full = sd(data$LeadersPreference) #1.097294
se_LeadersPreference_Full = sd(data$LeadersPreference)/sqrt(length(data$LeadersPreference)) #0.08067466

m1_LeadersPreference <- lm(LeadersPreference ~ PublicOpinion, data=data)
summary(m1_LeadersPreference) # p-value:  .683

## "The military schould carry out the will of the public"
PublicWillTable <-table(data$PublicWill)
prop.table(PublicWillTable)

Mean_PublicWill <- ddply(data, c("PublicOpinion"), summarise, mean = mean(PublicWill), sd = sd(PublicWill), sem = sd(PublicWill )/sqrt(length(PublicWill)))
print(Mean_PublicWill) #Public Supports: 2.06; Public Opposes: 2.35

Mean_PublicWill_Full <- (mean = mean(data$PublicWill)) #2.216216
sd_PublicWill_Full = sd(data$PublicWill) #1.091846
se_PublicWill_Full = sd(data$PublicWill)/sqrt(length(data$PublicWill)) #0.08027409

m1_PublicWill <- lm(PublicWill ~ PublicOpinion, data=data)
summary(m1_PublicWill) # p-value:  .0706

#----------------------
#Balance/Randomization Checks (Appendix, Section 2.1.1)

## Regress treatment condition variable on covariates (OLS). Used to produce Appendix, Table 4. 
BalanceCheck_PublicSupports<- lm(PublicSupports ~  Gender + Age + Mil_Branch + Combat + Status + Rank + Staff + Education + Race + Political_ID + Income, data=data) 
summary(BalanceCheck_PublicSupports)
BalanceCheck_PublicOpposes <- lm(PublicOpposes ~  Gender + Age + Mil_Branch + Combat + Status + Rank + Staff + Education + Race + Political_ID + Income, data=data) 
summary(BalanceCheck_PublicOpposes)

#Appendix, Table 4
stargazer(BalanceCheck_PublicSupports, BalanceCheck_PublicOpposes, title="OLS Models: Treatment Regressed on Covariates", align=TRUE)

#Balance Table of Covariate Means (Used to Produce Appendix, Table 3)

get_mean_sd <- function(vector){
  mean <- mean(vector, na.rm = T)
  sd <- sd(vector, na.rm = T)
  sem <- sd/(sqrt(length(vector)))
  output <- paste(round(mean, 3), " (", round(sem, 3), ")", collapse = "", sep = "")
  return(output)
}

balance_table <- group_by(data, PublicOpinion) %>%
  summarize(Branch = get_mean_sd(Mil_Branch),
            Rank = get_mean_sd(Rank),
            CombatExperience = get_mean_sd(Combat),
            StaffExperience = get_mean_sd(Staff),
            Gender = get_mean_sd(Gender),
            Race = get_mean_sd(Race),
            Age = get_mean_sd(Age),
            Education = get_mean_sd(Education),
            PoliticalID = get_mean_sd(Political_ID))


BalanceTable <- t(balance_table[,2:10])
colnames(BalanceTable) <- c("Public Supports Action", "Public Opposes Action")
stargazer(BalanceTable)

#----------------------
# Manipulation Check (Appendix Section 2.1.2, Table 5)

sum(data$PublicOpinion == 1 & data$Check == 2)/sum(data$PublicOpinion == 1) #0.9764706
sum(data$PublicOpinion == 2 & data$Check == 3)/sum(data$PublicOpinion == 2) #0.99
#------------------
#Demographics (Appendix, Section 1.3)

###Sample
##Gender 
GenderTable <-table(data$Gender)
prop.table(GenderTable)

##Rank
RankTable <-table(data$Rank)
prop.table(RankTable)

##Military Branch
BranchTable <-table(data$Mil_Branch)
prop.table(BranchTable)

##Service Headquarters/Joint Staff Experience 
StaffTable <-table(data$Staff)
prop.table(StaffTable)

##Combat Experience  
CombatTable <-table(data$Combat)
prop.table(CombatTable)

##Military Rank
CombatTable <-table(data$Combat)
prop.table(CombatTable)

##Education
EducationTable <-table(data$Education)
prop.table(EducationTable)

##Race
RaceTable <-table(data$Race)
prop.table(RaceTable)

##Political ID
PoliticalIDTable <-table(data$Political_ID)
prop.table(PoliticalIDTable)

###DOD Data (From DOD 2018 Demographics: Profile of the Military Community https://download.militaryonesource.mil/12038/MOS/Reports/2018-demographics-report.pdf)

# Percentage of Officers per service (Source: DOD 2018 Demographics, Page 13)--Note, the denominator of 212870 represents commissioned officers (see page 15, I exclude warrant officers)
#Army
92315/212870 #0.4336684

#Navy
54737/212870 #0.2571382

#Marine Corps
21332/212870 #0.09233673

#Air Force
62640/212870 #0.2711407

#Rank (Source: DOD 2018 Demographics, 15)

##O-1
28002/212870 #0.1315451

##0-2
28283/212870 #0.1328651

##0-3
74223/212870 #0.3486776

##0-4
43207/212870 #0.2029736

##0-5
27179/212870 #0.1276789

##0-6
11049/212870 #0.05190492

##0-7
440/212870 #0.002066989

# Gender Data (Source: DOD 2018 Demographics, 19)
# Race/Ethnicity (Source: DOD 2018 Demographics, 42)
# Education (Source: DOD 2018 Demographics, 42)
#---------------------------------------
# The additional analysis below was not included in the manuscript or appendix, but provides additional insights into the lack of heterogeneous effects 

# Mean Outcome (Recommendation on use of force) by covariates

##Branch Effects (Additional analysis not included in manuscript)
Branch_Table <- group_by(data, Mil_Branch) %>%
  summarize(Recommendation = get_mean_sd(Recommendation))

BranchTable <- t(Branch_Table[,2:2])
colnames(BranchTable) <- c("Army", "Navy", "Marines", "Air Force") #Recommendation Army: "1.159 (0.056)"; Navy: "1.059 (0.041)"; Marines: "1 (0)"; Air Force "1.097 (0.029)"

##Rank Effects (Additional analysis not included in manuscript)
Rank_Table <- group_by(data, Rank) %>%
  summarize(Recommendation = get_mean_sd(Recommendation))

RankTable <- t(Rank_Table[,2:2])
colnames(RankTable) <- c("O-1", "O-2", "O-3", "O-4", "O-5", "O-6", "O-6")

## Staff Effects (Additional analysis not included in manuscript)
meanRecommendation_Staff <- mean(data$Recommendation[data$Staff==2]) #1.071429
seRecommendation_Staff <- sd(data$Recommendation[data$Staff==2])/sqrt(length(data$Recommendation[data$Staff==2])) #0.0310041

meanRecommendation_NoStaff <- mean(data$Recommendation[data$Staff==1]) #1.121739
seRecommendation_NoStaff <- sd(data$Recommendation[data$Staff==1])/sqrt(length(data$Recommendation[data$Staff==1])) #0.03062488

Staff_Table <- group_by(data, Staff) %>%
  summarize(Recommendation = get_mean_sd(Recommendation))

StaffTable <- t(Staff_Table[,2:2])
colnames(StaffTable) <- c("No Staff Experience", "Staff Expereince")

## Combat Effects (Additional analysis not included in manuscript)
meanRecommendation_Combat <- mean(data$Recommendation[data$Combat==2]) #1.089286
seRecommendation_Combat <- sd(data$Recommendation[data$Combat==2])/sqrt(length(data$Recommendation[data$Combat==2])) #0.02706578

meanRecommendation_NoCombat <- mean(data$Recommendation[data$Combat==1]) #1.123288
seRecommendation_NoCombat <- sd(data$Recommendation[data$Combat==1])/sqrt(length(data$Recommendation[data$Combat==1])) #.03874558

## Gender Effects (Additional analysis not included in manuscript)
Gender_Table <- group_by(data, Gender) %>%
  summarize(Recommendation = get_mean_sd(Recommendation))

GenderTable <- t(Gender_Table[,2:2]) 
colnames(GenderTable) <- c("Male", "Female")

## Race Effects (Additional analysis not included in manuscript)
Race_Table <- group_by(data, Race) %>%
  summarize(Recommendation = get_mean_sd(Recommendation))

RaceTable <- t(Race_Table[,2:2])
colnames(RaceTable) <- c("White", "Black", "Asian", "Hawaiian/Pacific Islander", "Hispanic", "Mixed", "Other" )

## Education Effects (Additional analysis not included in manuscript)
Education_Table <- group_by(data, Education) %>%
  summarize(Recommendation = get_mean_sd(Recommendation))

EducationTable <- t(Education_Table[,2:2])
colnames(EducationTable) <- c("Bachelors", "Bachelors+", "Postgraduate Complete")

## Political ID Effects (Additional analysis not included in manuscript)
PoliticalID_Table <- group_by(data, Political_ID) %>%
  summarize(Recommendation = get_mean_sd(Recommendation))

PoliticalIDTable <- t(PoliticalID_Table[,2:2])
colnames(PoliticalIDTable) <- c("Very Liberal", "Liberal", "Moderate", "Conservative", "Very Conservative")
